% This script enables you to plot the torque vs RPM curve for a 
% variable speed generator defined by the parameters in Ashes
% If the default input values are kept, it will produce the curve
% for the NREL 5 MW wind turbine, as shown in figure 7-2 of the
% NREL 5+MW wind turbine specs (Jonkman 2009 Definition of a 5+MW 
% wind turbine for offshore development)

% Written by Loup, 2022-02-09

% Please do mention Simis AS if you use it.

clear all
close all


set(0, 'DefaultLineLineWidth', 1.5);


%% Input data, from Ashes
RPMr = 1173.7;        %rated RPM
Pr = 5000000;         %rated electrical power
eff = 0.944;          %generator efficiency

Kopt_perc = 0.8176;   %proportionality factor region 2
T115 = 0.570844;      %transition from region 1 to region 1.5
T152 = 0.74209;       %transition from region 1.5 to 2
T253 = 0.99;          %transition from region 2.5 to 3
slip = .1;            %generator slip


%% Computation of the different regions of the generator
Pmr = Pr/eff;                               %rated mechanical power
Tr = Pmr/(RPMr*pi/30);                      %rated torque

Kopt = Kopt_perc/(RPMr.^2/Tr);              %proportionality factor as defined in the NREL specs

RPM = [0:RPMr/1000:RPMr*1.2];                         %range of RPMs for the plot

RPM1 = [0,RPMr*T115];                       %range of RPMs of region 1
T1 = [0,0];                                 %range of torques of region 1

RPM15 = [RPMr*T115,RPMr*T152];              %range of RPMs of region 1.5
T15 = [0,Kopt*(RPMr*T152)^2];               %range of torques of region 1.5

a = (Tr/T253)/(T253*RPMr-(RPMr-slip*RPMr)); %slope of region 2.5
b = (Tr/T253) - a*T253*RPMr;                %y axis intercept of the curve of region 2.5


RPM225 = RPM(find(a*RPM+b>=RPM.^2*Kopt,1)); %RPM at transition between regions 2 and 2.5
T225 = RPM225.^2*Kopt;                      %torque at transition between regions 2 and 2.5



%% Plotting
plot(RPM,Kopt*RPM.^2,'k')
hold on
plot(RPM1,T1,'b')
plot(RPM15,T15,'b')
plot([RPM15(end):RPM225],[RPM15(end):RPM225].^2*Kopt,'b')
plot([RPM225 T253*RPMr],[T225 Tr/T253],'b')
plot([T253*RPMr:RPM(end)],30*Pmr./(pi.*[T253*RPMr:RPM(end)]),'b')

ylims = [0 RPM(end).^2*Kopt];
ylim(ylims)
xlim([0 RPM(end)])

plot([T115 T115]*RPMr,ylims,'k--','linewidth',0.5)
plot([T152 T152]*RPMr,ylims,'k--','linewidth',0.5)
plot([RPM225 RPM225],ylims,'k--','linewidth',0.5)
plot([T253 T253]*RPMr,ylims,'k--','linewidth',0.5)



legend({'Optimal','Variable speed controller'},'location','west')%,'fontsize',10)

xlabel('Generator speed [RPM]')
ylabel('Torque Nm')
grid on
set(gca,'FontSize',15)